* Nybom and Stuhler, Interpreting Trends in Intergenerational Mobility - Figure A4
* Simulate moments for a simple compulsory schooling reform

* Prepare Stata
clear
set more off 
set obs 1000
cap log close
log using MinimumSchooling_AnalyticMoments_Generations.log , replace


** (I) CHOOSE MODEL AND ITS PARAMETERS

* We choose parameters for the three structural equations we have in the manuscript, plus intercept, plus compulsory schooling parameter
* Outcomment line 92 below to simulate decrease in returns on the labor market

* Third structural equation (variance set to one)
local lambda = 0.6
local var_ue = 1-`lambda'^2
* second structural equation
local alpha_h = 0
local gamma_h = 1
local theta = 2 
local var_uh = 3
* first structural equation   
local alpha_y = 9
local delta = 0.2
local gamma_y = 0
local var_uy = 0.1
* minimum schooling limit
local complimit = 9 

* compute summarized parameter for 2-equations model
local gamma = `gamma_y' + `delta'*`gamma_h'
di "Gamma: " `gamma'
local rho = `delta'*`theta'
di "Rho: " `rho'
* compute expectations and variances in pre-shock steady state (solved on paper)
local exp_y = (`alpha_y'+`delta'*`alpha_h')/(1-`gamma_y'-`delta'*`gamma_h')
di "E[y]: " = `exp_y'
local exp_h = `alpha_h' + `gamma_h'*`exp_y'
di "E[h]: " = `exp_h'
local var_y = ((`rho'^2) + (2*`rho'^2*`gamma_h'*`lambda'*`delta'+2*`rho'^2*`gamma_y'*`lambda')/(1-`gamma'*`lambda') + `delta'^2*`var_uh' + `var_uy') / (1-`gamma'^2-`gamma_y'*`rho')
di "Var(y): " `var_y'					
local var_h = `gamma_h'^2*`var_y' + `theta'^2 + (2*`theta'^2*`gamma_h'*`lambda'*`delta')/(1-`gamma'*`lambda') + `var_uh'
di "Var(h): " `var_h'	
local cov_ye = `rho'/(1-`gamma'*`lambda')
di "Cov(y,e): " `cov_ye'
local beta_inc = `gamma' + ((`rho'^2*`lambda')/(1-`gamma'*`lambda'))/`var_y'
di "Beta_inc: " `beta_inc'
local beta_educ = (`gamma_h'^2*`beta_inc'*`var_y'+(`theta'*`gamma_h'*`lambda'^2+`theta'*`gamma_h')*`cov_ye'+`theta'^2*`lambda' + `delta'*`var_uh') / `var_h' 
di "Beta_educ: " `beta_educ'


** (II) SIMULATE DATA
set obs 2000000

* variables to be stored: assume shock hits in generation 16
forvalues i=15(1)19 {
	gen coeff_h`i'=.
	gen cov_h`i'=.
	gen corr_h`i'=.
	gen var_h`i'=.
	gen coeff_y`i'=.
	gen cov_y`i'=.
	gen corr_y`i'=.
	gen var_y`i'=.
}

forvalues i=1(1)5 {   // loop: simulate data and estimate coefficients multiple times

	** data generating process
	
	* generate/simulate first generation
	gen e1=rnormal()
	gen h1=rnormal()
	gen y1=rnormal()
	* generate pre-reform generations (verify that moments converge to pre-shock steady state)
	forvalues t=2(1)15 {	
		local lt = `t'-1
		qui gen e`t' = `lambda'*e`lt' + rnormal(0,sqrt(`var_ue'))
		qui gen h`t' = `alpha_h' + `gamma_h'*y`lt' + `theta'*e`t' + rnormal(0,sqrt(`var_uh'))
		qui gen y`t' = `alpha_y' + `gamma_y'*y`lt' + `delta'*h`t' + rnormal(0,sqrt(`var_uy'))
	} 
	* generate reform generations 
	forvalues t=16(1)20 {  	
		local lt = `t'-1
		qui gen e`t' = `lambda'*e`lt' + rnormal(0,sqrt(`var_ue'))
		qui gen h`t' = `alpha_h' + `gamma_h'*y`lt' + `theta'*e`t' + rnormal(0,sqrt(`var_uh'))
		qui gen h`t'cens = max(h`t',`complimit')
		qui gen y`t' = `alpha_y' + `gamma_y'*y`lt' + `delta'*h`t'cens + rnormal(0,sqrt(`var_uy'))	// COMMENT THIS LINE AND
		*qui gen y`t' = `alpha_y' + `gamma_y'*y`lt' + (`delta'-0.02)*h`t'cens + rnormal(0,sqrt(`var_uy'))	// OUTCOMMENT THIS LINE TO SIMULATE DECREASE IN RETURNS ON THE LABOR MARKET
	}

	* compute simluated moments
	gen h14cens  = h14  // to simplify code below
	gen h15cens  = h15  
	forvalues son=15(1)19 {
		local father = `son'-1
		
		* moments: h or h_censored
		reg h`son'cens h`father'cens , nohe
		replace coeff_h`son'=_b[h`father'cens] in `i'
		
		correlate h`son'cens h`father'cens , covariance
		mat A=r(C)
		replace cov_h`son'=A[2,1] in `i'
	
		correlate h`son'cens h`father'cens 
		mat A=r(C)
		replace corr_h`son'=A[2,1] in `i'
	
		su h`son'cens
		replace var_h`son'=r(Var) in `i'
	
		* moments: y
		reg y`son' y`father' , nohe
		replace coeff_y`son'=_b[y`father'] in `i'
	
		correlate y`son' y`father' , covariance
		mat A=r(C)
		replace cov_y`son'=A[2,1] in `i'
	
		correlate y`son' y`father' 
		mat A=r(C)
		replace corr_y`son'=A[2,1] in `i'
	
		su y`son'
		replace var_y`son'=r(Var) in `i'
	 	
	}
	drop e* h* y*
}

su coeff_h* coeff_y* corr_h* corr_y* cov_h* cov_y* var_h* var_y*


